Predicting the spin-lattice order of frustrated systems from first-principles 
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A novel general method of describing the spin-lattice interactions in magnetic solids was proposed 
in terms of first principles calculations. The spin exchange and Dzyaloshinskii-Moriya interactions 
cis well as their derivatives with respect to atomic displacements can be evaluated efficiently on the 
basis of density functional calculations for four ordered spin states. By taking into consideration 
the spin-spin interactions, the phonons, and the coupling between them, we show that the ground 
state structure of a representative spin-frustrated spinel, MgCr204, is tetragonally distorted, in 
agreement with experiments. However, our calculations find the lowest energy for the collinear spin 
ground state, in contrast to previously suggested non-coUinear models. 
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I. INTRODUCTION 

Spin frustrated systemsi""— have recently attracted con- 
siderable attention because of their novel magnetic prop- 
erties. The geometrically frustrated spin lattice gener- 
ally leads to numerous degenerate spin configurations. 
In principle, a strongly frustrated system mostly has 
no long-range spin order, but the degeneracy can of- 
ten be lifted by the spin-lattice coupling; a symmetry- 
lowering lattice distortion gives rise to a long-range mag- 
netic order at low temperature. In the emerging field of 
multiferroics^ it was found that many frustrated mag- 
nets (such as RMnOg^ and RMn20g^) display a large 
magnetoelectric coupling. Magnetic frustration com- 
bined with a striking spin-lattice coupling appears to 
cause their multiferroic properties. In particular, the ex- 
change striction (a kind of spin-lattice coupling arising 
from the dependence of the symmetric exchange interac- 
tions on the atomic positions) is considered to induce the 
ferroelectricity in some collinear antiferromagnets such as 
RMn205<i An alternative mechanism of multiferroicity 
comes from the inverse effect of Dzyaloshinskii-Moriya 
(DM) interaction^ (another kind of spin-lattice coupling 
arising from the dependence of antisymmetric exchange 
interactions on the atomic positions) , which explains the 
ferroelectricity in many noncoUinear spiral magnets such 
as RMnOg^"— and Ni3V208i--l A significant contribution 
of the symmetric Si • Sj-type magnetostriction to the fer- 
roelectricity in RMnOa was also revealed by Mochizuki 
et ai, who estimated the dependence of the nearest- 
neighbor (NN) ferromagnetic (FM) coupling in RMnOa 
on the Mn-O-Mn bond angle in an empirical way^iS, It 
is clear that a quantitative discription of the spin-lattice 
coupling is desirable for further study of frustrated mag- 
nets. 

In their pioneering work, Fennie and Rabe^ developed 
a first principles method to calculate the second order 
derivatives of the symmetric exchange parameter and the 



spin-phonon coupling parameter, and they studied the 
influence of magnetic order on the optical phonons of 
the geometrically frustrated spinel ZnCr204. However, 
their work does not fully describe how the spins cou- 
ple to the lattice^i^, and how to predict the spin-lattice 
order in frustrated systems from first principles is un- 
clear. Here in this work, we propose a general first princi- 
ples description of spin-lattice coupling, which allows one 
to efficiently evaluate the symmetric and antisymmet- 
ric exchange interaction parameters and their first-order 
derivatives with respect to atom displacements. As a test 
of our method, we examined a representative frustrated 
system, the spinel MgCr204, to find that its ground-state 
structure is tetragonally distorted with collinear antifer- 
romagnetic (AFM) spin configuration. 

In section II, we will describe a method for computing 
the spin-lattice coupling parameters. The approach for 
predicting the spin-lattice ground state using the spin- 
lattice coupling parameters will be discussed in section 
HI. Then we will apply the methods to find the spin- 
lattice ground state of MgCr204 in section IV. Finally, 
we will summarize our work in section V. 



II. QUANTITATIVE DESCRIPTION OF THE 
SPIN-LATTICE COUPLING 

^. Method for computing the symmetric exchange 
parameters 



Consider a classical Heisenberg spin system, whose en- 
ergy can be written as E = Eq -\- Espt 



Din 1 where E^p^n 



Sj is the spin exchange term with |Si| =5*, 



and Eq is the energy due to other interactions (e.g., the 
lattice elastic energy), which depends on the atom dis- 
placements but not the spin orientations. The spin ex- 
change interactions are short range interactions and be- 
come negligible when the distance between the spin sites 
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i and j is longer than a certain critical distance Re- Given 
a supercell large enough so that any spin site has no inter- 
action with its neighboring cells, we extract the exchange 
interactions as follows. Without loss of generality, con- 
sider a particular exchange interaction J12 between spin 
sites 1 and 2. The spin Hamiltonian can be written as: 

Espin — Jl2^1 • S2 + Si • Ki + S2 ■ K2 + Bother: (1) 
where Ki — J2i^l.2'^'^i^i' = 'l2i^l.2 "^'^^^i^ Eother = 

Si j=£i 2 Jij^i ■ ^j- It should be noted that Ki, K2, and 
Eother do not depend on the spin directions of sites 1 and 
2. Consider the following four collinear spin states (with 
z as the spin quantization axis): (1) Sf — 8,82— S; (2) 
S! = S, Si = -S; (3) Sf = -S, Si = S; (4) Sf = -S, 
Si = —S. In these four spin states, the spin orientations 
for the spin sites other than 1 and 2 are the same. One 
can easily show that the four states have the following 
energy expressions: 



Ei=Eo+ Bother + Jl2S^ + KiS + K2S, 
E2=Eo+ Eother - Jl2^' + K^S - K2S, 
E3=Eo+ Eother " Jl2^' - K^S + K2S, 
E4 = Eq + Eother + Jl2S^ ~ KiS ~ K2S. 



Thus, J12 is extracted by the formula: 



J12 



El + E4 ~ E2 



E. 



452 



(2) 



(3) 



The total energies of the four states can be calculated 
using density functional theory (DFT). The method de- 
scribed above is a kind of mapping analysis^^ that has 
been used widely to extract the exchange parameters. In 
the usual mapping process, one usually considers the spin 
states with small supercells to reduce computational de- 
mand. Here the use of a large supercell has the advantage 
that the extraction of a particular exchange interaction is 
independent of other exchange interactions and all total 
energies are computed using the same supercell. In prin- 
ciples, the accuracy of the exchange parameters from this 
method is only limited by the predefined cutoff distance 
Rc, which can be checked by increasing the supercell size. 



B. Method for computing the derivatives of the 
symmetric exchange parameters 

The above approach considered how to evaluate the 
exchange parameters for a given structure. We now 
examine the dependence of the total energy on the 
atom displacements: i5(ui, u„. Si, Sm) = i?o(U) + 
Y.<2,j> ■ where U = (ui,...,u„) and Ufc 

(1 < A: < n) denote the displacements of atom k from 
a reference structure, n and m are the total number of 
atoms and total number of spin sites in the supercell, re- 
spectively. By taking the partial derivative of the above 
equation with respect to Uka, we obtain: 



dE dEo 



duh 



duk 



E 

<i,3> 



dJi, 



du 



■Si ■ Sj 



koL 



(4) 



In terms of the same four spin states used for extracting 
J12, the derivative of J12 with respect to Uka is found as: 



dJi2 _ 1 , dEi dEi 



dE2 dEs . 

dUka iS'^^dUka ' duka duka duka' 



(5) 



Here 



dEj 
duka 



{i = 1, ...,4) is the force acting on the atom 



k along the a direction. The force can be computed us- 
ing the Hellmann-Feynman theorem and is readily avail- 
able in many standard DFT schemes. From Eq. [SJ we 
can see that the dependence of the exchange parameter 
J12 on all the atom displacements can be computed by 
performing four static total energy calculations. This 
means that the calculation of first order derivative of 
the exchange parameter does not require extra calcula- 
tions if one calculates the exchange parameter using our 
above method. Therefore, our new approach utilizing the 
Hellmann-Feynman forces has a great computational ad- 
vantage over the finite difference method in which each 
requires several total energy calculations. 



C. Methods for computing DM interaction 
parameters, single-ion anisotropy parameters, and 
their derivatives 

Our method for calculating the symmetric exchange 
parameter and its derivative can be also extended to 
compute the antisymmetric DM interaction parameter 
( D) and single-ion anisotropy (SIA) parameter {A) and 
their derivatives. Let us describe the method of calcu- 
lating the DM interaction parameter (vector D12) be- 
tween spin site 1 and spin site 2. Here, we calculate 
the three components Df2 , 0^2 ^^^d ^12 of the DM vec- 
tor separately for a general system, although sometimes 
the direction of the DM vector can be determined by 
the crystal symmetry. Without loss of generality, let us 
focus on the calculation of £>f2- We consider the fol- 
lowing four spin configurations in which the spins 1 and 
2 are oriented along the x— and y— axes, respectively: 
(1) Si = (5,0,0), S2 = (0,5,0), (2) Si = (5,0,0), 
S2 - (0,-5,0), (3) Si = (-5,0,0), S2 = (0,5,0), (4) 
Si = (-5,0,0), S2 = (0,-5,0). In these four spin con- 
figurations, the spins of all the other spin sites are the 
same and are along the z direction. The spin interaction 
energy for the four spin configurations can be written as: 



Espin — -Di25i 5|— 5i £'j'j5f-f5| Di^S^+Eg-^ 
As in the case of the symmetric exchange, we have 



ther ■ 



(6) 



DI2 = j^{Ei+Ei~E2-E3), 

ODI2 ^ 1 ( dEi I dEj _ dE2 _ dE3 \ 



(7) 



Since the DM interaction is a consequence of spin-orbit 
coupling (SOC), it is necessary that the energies of the 
four ordered spin states be determined by DFT calcu- 
lations with SOC effects taken into consideration. The 
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other two components and can be computed in 
a similar manner. This is a general method to compute 
both the direction and the magnitude of a DM inter- 
action parameter. For the calculation of the single-ion 
anisotropy parameter, we consider the spin site 1. If 
the spin has an easy-axis (with local z' axis) or easy- 
plane anisotropy, the single-ion anisotropy term can be 
expressed as: iJsm = AiS^,. To evaluate Ai, we con- 
sider the four spin states in which the spin directions 
for site 1 are along z', —z', x\ and —x' with the spins 
at all the other spin sites along the y direction. One 
can easily find that Ai = E,+E2-E3-Ei ^ 
1 idE^ ^ dE^ _ jm^ _ dEA.), Here the totaT en- 



ergy calculations should include SOC effects because the 
single-ion anisotropy is a consequence of SOC. 



III. PREDICTION OF THE SPIN-LATTICE 
ORDER 



It was shown that the spin- lattice coupling may lead to 
a distortion of the lattice to lower the exchange energy 
and relieve the frustration of a frustrated spin system 
(the "spin- Teller" effecti^). How the lattice distorts is 
determined by the balance between the spin exchange 
and lattice elastic energies, and depends on the param- 
eters of the full spin-lattice coupled Hamiltonian. With 
the exchange parameters and their first order derivatives 
with respect to the atom displacements in hand, one can 
predict its spin-lattice ground state with a lower crystal 
symmetry. At high temperature, a magnetically frus- 
trated system is usually in a disordered paramagnetic 
(PM) state with a high symmetry. We now write the to- 
tal energy of a frustrated system with the PM state as a 
reference state: 



-E'(ui, u„, Si, 



E 



PM 



Eph + Eg 



(8) 



where 



ph 



1/2 E 



,Uj(j is the phonon Hamil- 



tonian (C"^ is the force constant), and Eg. 



E<..,>[^.{u)s,.s, 



Here Ji 



-D,;,(U)-(S,xS,)] + E,A,(U)Sf 



7° 



9Ji 



computed using the PM structure. We have similar ex- 
pressions for Dy (U) and Ai{XJ). The particular lattice 
distortion leading to the lowest energy for a given spin 
configuration can be obtained by solving the following 
linear equations: 



7° 



and 



iz' ■ 

are 



dE 
du 



ka 



E 
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Ck3 '^jP 



dE 



spin 



duka 



= 



(9) 



where 1 < k < n and a — 1,2,3. The lowest energy for a 
given spin configuration can then be calculated by using 
Eq. [51 The spin-lattice ground state can be found by 
comparing the energies of different spin configurations. 




FIG. 1: (Color online) (a) Structure of MgCr204 and ex- 
change paths, (b) and (c) show the side and top views of the 
derivative of the NN exchange Ji with respect to the atom 
displacements. 



IV. APPLICATION TO MGCR2O4 
A. Computational details of DFT calculations 

Our total energy calculations are based on the DFT 
plus the on-site repulsion (U) method^^ within the local 
density approximation (LDA-I-U) on the basis of the pro- 
jector augmented wave method^^ encoded in the Vienna 
ab initio simulation package-?. We used the on-site repul- 
sion [/ = 3 eV and the exchange parameter J = 0.9 eV 
on Cr, which reproduce the dominant features of the pho- 
toemission and band gap data in sulfur Cr^"'" spinelsi^. 
The plane-wave cutoff energy was set to 400 eV. 



B. Exchange parameters and their derivatives 

We now apply the above method to calculate the ex- 
change parameters and their derivatives of MgCr204^- 
to find its spin-lattice coupled ground state. MgCr204, 
which crystallizes in a normal spinel structure with the 
Cr-^^ ions {S = 3/2) forming a pyrochlore lattice [see 
Fig-UJa)], is strongly frustrated due to the strong AFM 
interactions between the NN Cr'^+ spins. We consider 
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FIG. 2: (Color online) The ground state of MgCr2 04 by con- 
sidering the full spin-lattice coupled Hamiltonian. The num- 
bers (in A) denote the Cr-Cr distances. The exchange pa- 
rameter for the Cr-Cr pair with the 2.911 A distance is 5.47 
meV, while the exchange parameter for the Cr-Cr pair with 
the 2.937 A distance is 3.84 meV. 



first the MgCr204 structure optimized with the FM spin 
state. The optimized lattice constant is 8.277 A. We 
calculate all symmetric spin exchange interactions up to 
the fourth NN pairs [see Fig. [Ua)] using a 2 x 2 x 2 su- 
percell of the MgCr204 conventional cubic cell; the NN 
exchange Ji within each Cr4 tetrahedron and the farther 
NN exchanges J2, J3, J4, and J5. We find that Ji =4.56 
meV, J2 = 0.01 meV, J3 = 0.26 meV, J4 = 0.08 meV, 
and J5 = —0.01 meV. The NN exchange Ji is strongly 
AFM while all next NN exchange interactions are almost 
negligible. 

Using Eq. [5l we calculate the derivatives of the ex- 
change parameters with respect to atom displacements 
using the optimized FM structure. Our results show the 
NN exchange Ji depends strongly on the positions of 
those atoms shown in Figs. [IJb) and (c). In particular, 
the largest derivative of the exchange interaction 

between two NN Cr ions occurs when site k is one of 
two Cr ions (We hereafter refer to this vector as Jic,.)- 
When the two Cr ions move close to each other, the NN 
exchange Ji increases. The magnitude of the derivative 
is as large as 43.40 meV/A, which is close to the value 
(40.25 meV/A) extracted by the finite difference method. 
This can be understood because when the distance be- 
tween two Cr ions becomes short, the direct overlap be- 
tween the t2g orbitals of the two Cr^+ {(P) ions becomes 
stronger. The NN exchange Ji also depends substantially 
on the positions of the two bridging oxygen atoms with 
its derivative approximately along the direction from the 
midpoint of the two NN Cr ions to the oxygen ion (here- 
after this vector is referred to as J^^q)- We find the direc- 
tion of is due to the fact that the anti-bonding re- 
pulsion between O p orbital and Cr t2g orbitals becomes 
weaker and the t2g orbitals of the two Cr^+ (cfi) ions can 
have a better overlap if the bridge O atoms move away 
from the Cr-Cr pair. Another possible explanation is that 
the increased Cr-0 distance might result in a smaller FM 
contribution, enhancing the overall AFM coupling. The 
derivatives of other symmetric exchange parameters are 
found to be vanishingly small. 

It was suggested that the DM interaction might be 
relevant to the spin-lattice order in a similar system.^" 



As expected from the symmetry analysis, the DM vector 
for a Cr-Cr edge of each Cr4 tetrahedron is perpendicular 
to the Cr-Cr bond and is parallel to the opposite edge 
of the Cr4 tetrahedron. Using our method, we find the 
magnitude of the DM parameter to be 0.03 meV (0.7% of 
Ji). We have checked that our method is accurate enough 
to predict reliable DM vectors. And the derivative of the 
DM parameter with respect to the Cr ion position of the 
Cr-Cr pair has the largest magnitude of 0.41 me V/A. Our 
calculations show that the Cr^+ ion has an easy plane 
anisotropy with the plane perpendicular to the three-fold 
rotational axis z' . The calculated SIA parameter is —0.05 
meV (1.1% of Ji) and the largest derivative of the SIA is 
0.18 meV/A. Our first principles calculation established 
that, for MgCr204, the DM parameter, SIA parameter, 
and theirs derivatives are negligible compared to NN Ji . 



C. Spin- lattice ground state of MgCr2 04 

At high temperature, MgCr204 is paramagnetic (FM) 
with a cubic symmetry. To simulate the PM state, we 
use a special quasirandom structure (SQS)^^ spin config- 
uration. The magnetic unit cell of the SQS spin config- 
uration is four times as large as the chemical unit cell of 
MgCr204. Then we fix the lattice constant to that (8.259 
A) of the SQS structure and relax the internal atomic co- 
ordinates with the FM spin configuration. In this way, 
we can get approximate exchange parameters J°- (4.80 

meV), their derivative (|JiCil ~ 49.11 meV/A and 

|J'^q| — 25.15 meV/A), and the force constants C"j^ ap- 
propriate for the PM state. The differences between the 
force constants of the FM state and those of the PM 
state are small because they are of second order effect, 
and thus are neglected in this work. All these parameters 
are necessary in finding the spin-lattice coupled ground 
state of MgCr204 by solving Eq. [9] The force constants 
are obtained by finite difference as in phonon calculations 
by the direct method. The force constants are considered 
fixed and independent of the atomic displacement. 

Our calculations show that the dominant exchange in- 
teraction in MgCr204 is the NN AFM symmetric ex- 
change interaction Ji. For the Heisenberg Hamiltonian 
with the NN exchange interaction Ji on the pyrochlore 
lattice {Hnn = Z]<NNj,j> ^^i^^ ' ^i)' ^^'^ degeneracy of 
the spin ground state is macroscopic: If for any of the 
Cr tetrahedra, the sum of the four spins is zero, then it 
is a spin ground statcji^ We will use our calculated pa- 
rameters (not only the exchange parameter, but also its 
derivatives) and solve the full spin-lattice coupled Hamil- 
tonian to determine the ground state of MgCr204. As 
a first step, we consider the case where the spin-lattice 
ground state has the same size as the primitive chemi- 
cal unit cell. In this case, we can easily generate a spin 
configuration^^ that is one of the highly degenerate spin 
ground states of the Heisenberg Hamiltonian H^n. We 
generate two thousand of such spin configurations and 
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calculate the energy of the spin-lattice coupled system, 
to find that the spin configuration shown in Fig.[2]has the 
lowest energy and thus is the ground state of the spin- 
lattice coupled Hamiltonian. The spin state is collinear 
with two up spins and two down spins. To confirm the 
above prediction from the model Hamiltonian analysis, 
we carry out DFT calculations to optimize both the lat- 
tice parameters and the internal coordinates of MgCr204 
with the collinear AFM spin state. The relaxed (i.e., 
distorted) structure is calculated to have a lower energy 
by 6.33 meV/Cr than the unrelaxed structure with the 
same collinear AFM spin state. In the relaxed struc- 
ture, the distance between NN spin up Cr'^^ ion and spin 
down Cr'^^ ion is smaller by 0.026 A than that between 
NN Cr'^^ ions with the same spin direction. The ex- 
change parameters for the two different Cr-Cr exchange 
interactions are 5.47 meV and 3.84 meV, respectively, 
to be compared with the value (4.80 meV) for the undis- 
torted structure. This difference is due to the spin-lattice 
coupling, which makes the exchange parameter between 
AFM (FM) coupled spins larger (smaller). The relaxed 
MgCr204 structure is tetragonal (space group I^i/amd, 
No. 141) with a = 5.873 A and c 8.160 A. Exper- 
imentally, the spinel MgCr204 was found to undergo a 
sharp first order transition at ~ 12.4 K from a cu- 
bic paramagnetic phase (space group FdZm) to a tetrag- 
onal antiferromagnetically ordered structure {lAi/amd, 
a = 5.8961A and c = 8.3211 A at 10 K)}^ Our first prin- 
ciples result thus confirms the lAi/amd space group of 
the ground state of MgCr204 and c < V^a below T^- It 



should be noted that our work predicts a collinear AFM 
ground state with the propagation vector q — (0, 0, 0) 
with respect to the tetragonal lattice, in contrast to the 
previously proposed non-collinear magnetic modelsJ^i^^ 
This calls for further ultra-low temperature experiments 
to verify our prediction. 



V. CONCLUSION 

In summary, we proposed a general and efficient 
method to quantitatively describe spin-lattice coupling. 
This method allows one to evaluate each spin exchange 
(and DM interaction) as well as its derivatives with re- 
spect to atom displacements on the basis of DFT cal- 
culations for four ordered spin states. By applying 
this method to the spin- frustrated spinel MgCr204, we 
showed that it undergoes a structural transition from 
the cubic to a tetragonal structure with collinear AFM 
spin configuration. Our method provides an efficient first 
principles way of describing the interplay between spin 
order and lattice distortion in frustrated magnetic sys- 
tems. 
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